
*** Replication do file for Linking Local Infra. paper
*** By: Christian Baehr, Ariel BenYishay, and Bradley Parks
*** Published in the Journal of the Association of Environmental and Resources Economists
*** Created October 23, 2020

* ssc install reghdfe
* ssc install outreg2

* set your working directory
cd "/Users/christianbaehr/Downloads/cambodia_replication"

* check that required datasets are in working directory
local files "main_panel pid commune_panel"
foreach i in `files' {
	confirm file "`i'.dta"
}

********** Load and Process Data **********

* load the panel dataset
use main_panel, clear

* set up the panel with "cell_id" as the panel variable and "year" as the time variable
xtset cell_id year

* aggregate the road project count. The aggregated measure will give the total number of completed road projects within 5km of grid cell x
egen trt_overall_road = rowtotal(trt_road_?kmband)
egen max_trt_overall_road = max(trt_overall_road), by(cell_id)

* aggregate the irrigation project count
egen trt_overall_irrigation = rowtotal(trt_irrigation_?kmband)
egen max_trt_overall_irrigation = max(trt_overall_irrigation), by(cell_id)
* aggregate the total project count (all project types)
egen trt_overall = rowtotal(trt_?kmband)

* create an "other" project count. This is just the total project count minus the number of road and irrigation projects
gen trt_overall_else = trt_overall - (trt_overall_road + trt_overall_irrigation)

* generate a variable with the maximum number of projects a cell is treated by
egen max_trt_overall = max(trt_overall), by(cell_id)
* drop cells which are never treated
drop if max_trt_overall==0
* generate a new categorical variable for max_treatment. Values above 20 are topcoded
egen cut_max_trt_overall = cut(max_trt_overall), at(0(1)20 1000) label

* generate a variable indicating whether a cell has received any treatments yet
gen trt_overall_pos = (trt_overall>0)
* identify the year of first treatment
gen yr_tmp = year if trt_overall_pos==1 & l1.trt_overall_pos==0
* populate all obs of each cell with the year of first treatment
egen year_first_pos = max(yr_tmp), by(cell_id)
* drop temp variable
drop yr_tmp

* generate a high population dummy
gen high_pop = (pop_density_2000>=1000)

* generate a variable indicating whether a cell has received any road-related treatments yet
gen trt_overall_road_pos = (trt_overall_road>0)
* identify the year of first road treatment
gen yr_tmp = year if trt_overall_road_pos==1 & l1.trt_overall_road_pos==0
* populate all obs of each cell with the year of first road-related treatment
egen year_first_road_pos = max(yr_tmp), by(cell_id)
* generate a "years to first road project variable". This idenfities the temporal distance of an observation from initial treatment
gen years_to_first_road = year - year_first_road_pos
* generate years since first treatment. All pre-treatment observations are zero
gen years_since_first_road = max(0, years_to_first_road) 
* drop temp variable
drop yr_tmp
* generate a new categorical variable for road treatment count. Values above 20 are topcoded
egen cut_trt_overall_road = cut(trt_overall_road), at(0(1)20 1000) label

* generate a variable indicating whether a cell has received any irrigation-related treatments yet
gen trt_overall_irrigation_pos = (trt_overall_irrigation>0)
* identify the year of first irrigation treatment
gen yr_tmp = year if trt_overall_irrigation_pos==1 & l1.trt_overall_irrigation_pos==0
* populate all obs of each cell with the year of first irrigation-related treatment
egen year_first_irrigation_pos = max(yr_tmp), by(cell_id)
* generate a "years to first irrigation project variable". This idenfities the temporal distance of an observation from initial irrigation treatment
gen years_to_first_irrigation = year - year_first_irrigation_pos
* generate years since first irrigation treatment. All pre-treatment observations are zero
gen years_since_first_irrigation = max(0, years_to_first_irrigation)
* drop temp variable
drop yr_tmp
* generate a new categorical variable for irrigation treatment count. Values above 20 are topcoded
egen cut_trt_overall_irrigation = cut(trt_overall_irrigation), at(0(1)20 1000) label

* generate indicator of whether a cell has received any non road or irrigation related treatments by year t
gen trt_overall_else_pos = (trt_overall_else>0)
* identify the year of first non road or irrigation treatment
gen yr_tmp = year if trt_overall_else_pos==1 & l1.trt_overall_else_pos==0
* populate all obs of each cell with the year of first non road or irrigation-related treatment
egen year_first_else_pos = max(yr_tmp), by(cell_id)
* generate a "years to first non road or irrigation project variable". This idenfities the temporal distance of an observation from initial non road or irrigation treatment
gen years_to_first_else = year - year_first_else_pos
* generate years since first non road or irrigation treatment. All pre-treatment observations are zero
gen years_since_first_else = max(0, years_to_first_else)

* compute the rowwise mean of yearly percent Seila funding variable across all years (1996-2003)
egen seila_total = rowmean(seila_pct_*)


********** Figures & Tables **********

*** Figure 1 

* save the panel in memory
preserve

* use the PID dataset for the first figure
use pid, clear

* histograms of treatment activity, disaggregated to rural transportation projects, irrigation projects, and other projects
twoway (histogram end_year if activity_type=="Rural Transport", discrete freq color(red%30) recast(scatter) start(2003) xlabel(2003(2)2019)) (histogram end_year if activity_type=="Irrigation", discrete freq color(blue%30) recast(scatter) start(2003) xlabel(2003(2)2019)) (histogram end_year if activity_type!="Rural Transport" & activity_type!="Irrigation", discrete freq color(green%30) recast(scatter) start(2003) xlabel(2003(2)2019)), xtitle("Year") ytitle("Projects completed") legend(order(1 "Roads" 2 "Irrigation" 3 "Other") rows(1)) bgcolor(white) graphregion(color(white)) saving(figure1, replace)

* restore the panel
restore


*** Table 1

* write out Table 1 summary statistics
outreg2 using table1.doc, replace sum(log) keep(treecover ndvi trt_1kmband trt_2kmband trt_3kmband trt_4kmband trt_5kmband temperature precipitation plantation_dummy concession_dummy protected_area_dummy ntl infant_mort distance_to_city distance_to_road distance_to_minecasualty pop_density_2000 bombing_dummy burial_dummy memorial_dummy prison_dummy)


*** Figure 3

* generate year labels for figure axes
loc year_lab ""
forv y = 2003/2018 {
	loc year_lab "`year_lab' `y'.year_first_pos#c.year = `y'"
}
* generate # of treatment labels for figure axes
forv t = 1/20 {
	loc t_lab "`t_lab' `t'.cut_max_trt_overall#c.year = `t'"
}

* a)
reghdfe ndvi ibn.year_first_pos#c.year if year<2003 [aw=cell_count_30m], absorb(cell_id) cluster(commune_id year)
loc ref = _b[2003.year_first_pos#c.year]
coefplot, keep(*.year_first_pos*) vertical yline(`ref') graphregion(color(white)) legend(off) xtitle("Year of first project") ytitle("Pre-2003 trend in NDVI") rename(`year_lab') xlabel(, labsize(vsmall) alternate) ylabel(-0.05(0.05)0.15) yscale(range(-0.05(0.05)0.15)) saving(figure3_a, replace)

* b)
reghdfe ndvi ibn.cut_max_trt_overall#c.year if year<2003 & cut_max_trt_overall>0 [aw=cell_count_30m], absorb(cell_id) cluster(commune_id year)
loc ref = _b[1.cut_max_trt_overall#c.year]
coefplot, keep(*.cut_max_trt_overall*) vertical yline(`ref') graphregion(color(white)) legend(off) xtitle("Total projects, 2003-2018") ytitle("Pre-2003 trend in NDVI") xlabel(1(1)20) ylabel(-0.05(0.05)0.15) yscale(range(-0.05(0.05)0.15)) saving(figure3_b, replace)

* c)
reghdfe treecover ibn.year_first_pos#c.year if year<2003 [aw=cell_count_30m], absorb(cell_id) cluster(commune_id year)
loc ref = _b[2003.year_first_pos#c.year]
coefplot, keep(*.year_first_pos*) vertical yline(`ref') graphregion(color(white)) legend(off) xtitle("Year of first project") ytitle("Pre-2003 trend in tree cover") rename(`year_lab') xlabel(, labsize(vsmall) alternate) ylabel(-0.02(0.01)0.01) yscale(range(-0.02(0.01)0.01)) saving(figure3_c, replace)

* d)
reghdfe treecover ibn.cut_max_trt_overall#c.year if year<2003 & cut_max_trt_overall>0 [aw=cell_count_30m], absorb(cell_id) cluster(commune_id year)
loc ref = _b[1.cut_max_trt_overall#c.year]
coefplot, keep(*.cut_max_trt_overall*) vertical yline(`ref') graphregion(color(white)) legend(off) xtitle("Total projects, 2003-2018") ytitle("Pre-2003 trend in tree cover") xlabel(1(1)20) ylabel(-0.02(0.01)0.01) yscale(range(-0.02(0.01)0.01)) saving(figure3_d, replace)


*** Table 2

* (1)
reghdfe ndvi trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, replace noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (2)
reghdfe treecover trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (3)
reghdfe ndvi i.high_pop#c.(trt_overall_road trt_overall_irrigation trt_overall_else) temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (4)
reghdfe treecover i.high_pop#c.(trt_overall_road trt_overall_irrigation trt_overall_else) temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (5)
reghdfe ndvi trt_overall_road years_since_first_road trt_overall_irrigation years_since_first_irrigation trt_overall_else years_since_first_else temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (6)
reghdfe treecover trt_overall_road years_since_first_road trt_overall_irrigation years_since_first_irrigation trt_overall_else years_since_first_else temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table2.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)


*** Table 3

* (1)
reghdfe ntl trt_overall_road trt_overall_irrigation trt_overall_else [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, replace noni nocons ctitle (Nighttime Light) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", N)

* (2)
reghdfe infant_mort trt_overall_road trt_overall_irrigation trt_overall_else [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, append noni nocons ctitle (Infant Mort.) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", N)

* (3)
reghdfe ndvi c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.(plantation concession protected_area) temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, append noni nocons ctitle (NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (4)
reghdfe treecover c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.(plantation concession protected_area) temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, append noni nocons ctitle (Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (5)
reghdfe ndvi c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.distance_to_minecasualty temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, append noni nocons ctitle (NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (6)
reghdfe treecover c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.distance_to_minecasualty temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table3.doc, append noni nocons ctitle (Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)


*** Figure 4

* quintiles of 30m cell counts. 30m cell count is the number of 30m grid cells involved in the aggregation to the 1km cell.
xtile q_count = cell_count_30m, nq(5)

reghdfe ndvi ibn.q_count#c.(trt_overall_road trt_overall_irrigation trt_overall_else) temperature precipitation, absorb(cell_id year) cluster(commune_id year)

* a)
coefplot, keep(*.q_count#c.trt_overall_road) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Quintile of share initially forested") ytitle("Effect of roads on NDVI") rename(1.q_count#c.trt_overall_road = 1 2.q_count#c.trt_overall_road = 2 3.q_count#c.trt_overall_road = 3 4.q_count#c.trt_overall_road = 4 5.q_count#c.trt_overall_road = 5) ylabel(-0.002(0.002)0.008) yscale(range(-0.002(0.002)0.008)) saving(figure4_a, replace)

* b)
coefplot, keep(*.q_count#c.trt_overall_irrigation) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Quintile of share initially forested") ytitle("Effect of irrigation on NDVI") rename(1.q_count#c.trt_overall_irrigation = 1 2.q_count#c.trt_overall_irrigation = 2 3.q_count#c.trt_overall_irrigation = 3 4.q_count#c.trt_overall_irrigation = 4 5.q_count#c.trt_overall_irrigation = 5) ylabel(-0.002(0.002)0.008) yscale(range(-0.002(0.002)0.008)) saving(figure4_b, replace)

reghdfe treecover ibn.q_count#c.(trt_overall_road trt_overall_irrigation trt_overall_else) temperature precipitation, absorb(cell_id year) cluster(commune_id year)

* c)
coefplot, keep(*.q_count#c.trt_overall_road) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Quintile of share initially forested") ytitle("Effect of roads on tree cover") rename(1.q_count#c.trt_overall_road = 1 2.q_count#c.trt_overall_road = 2 3.q_count#c.trt_overall_road = 3 4.q_count#c.trt_overall_road = 4 5.q_count#c.trt_overall_road = 5) ylabel(-0.01(0.01)0.02) yscale(range(-0.01(0.01)0.02)) saving(figure4_c, replace)

* d)
coefplot, keep(*.q_count#c.trt_overall_irrigation) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Quintile of share initially forested") ytitle("Effect of irrigation on tree cover") rename(1.q_count#c.trt_overall_irrigation = 1 2.q_count#c.trt_overall_irrigation = 2 3.q_count#c.trt_overall_irrigation = 3 4.q_count#c.trt_overall_irrigation = 4 5.q_count#c.trt_overall_irrigation = 5) ylabel(-0.01(0.01)0.02) yscale(range(-0.01(0.01)0.02)) saving(figure4_d, replace)


*** Figure 5

* a)
reghdfe ndvi trt_road_?kmband temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
coefplot, keep(trt_road_?kmband) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Distance band (km)") ytitle("Effect of roads on NDVI") rename(trt_road_1kmband = 0-1 trt_road_2kmband = 1-2  trt_road_3kmband = 2-3  trt_road_4kmband = 3-4  trt_road_5kmband = 4-5) ylabel(-0.005(0.005)0.01) yscale(range(-0.005(0.005)0.01)) saving(figure5_a, replace)

* b)
reghdfe ndvi trt_irrigation_?kmband temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
coefplot, keep(trt_irrigation_?kmband) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Distance band (km)") ytitle("Effect of irrigation on NDVI") rename(trt_irrigation_1kmband = 0-1 trt_irrigation_2kmband = 1-2  trt_irrigation_3kmband = 2-3  trt_irrigation_4kmband = 3-4  trt_irrigation_5kmband = 4-5) ylabel(-0.005(0.005)0.01) yscale(range(-0.005(0.005)0.01)) saving(figure5_b, replace)

* c)
reghdfe treecover trt_road_?kmband temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
coefplot, keep(trt_road_?kmband) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Distance band (km)") ytitle("Effect of roads on tree cover") rename(trt_road_1kmband = 0-1 trt_road_2kmband = 1-2  trt_road_3kmband = 2-3  trt_road_4kmband = 3-4  trt_road_5kmband = 4-5) ylabel(0(0.005)0.02) yscale(range(0(0.005)0.02)) saving(figure5_c, replace)

* d)
reghdfe treecover trt_irrigation_?kmband temperature precipitation [aw=cell_count_30m], absorb(cell_id year) cluster(commune_id year)
coefplot, keep(trt_irrigation_?kmband) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Distance band (km)") ytitle("Effect of irrigation on tree cover") rename(trt_irrigation_1kmband = 0-1 trt_irrigation_2kmband = 1-2  trt_irrigation_3kmband = 2-3  trt_irrigation_4kmband = 3-4  trt_irrigation_5kmband = 4-5) ylabel(-0.005(0.005)0.01) yscale(range(-0.005(0.005)0.01)) saving(figure5_d, replace)


*** Table 5

* preserve the 1km panel
preserve

* generate year dummies
cap qui tab year, gen(yr_dum)

* create copies of treatment variables for IV estimation
clonevar tor = trt_overall_road
clonevar toi = trt_overall_irrigation
clonevar toe = trt_overall_else

* (1)
ivreg2h ndvi temperature precipitation yr_dum* (tor toi toe=) [aw = cell_count_30m], fe cluster(commune_id) 
outreg2 using table5.doc, replace drop(yr_dum*) noni nocons ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (2)
ivreg2h treecover temperature precipitation yr_dum* (tor toi toe=) [aw = cell_count_30m], fe cluster(commune_id) 
outreg2 using table5.doc, append drop(yr_dum*) noni nocons ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* restore the 1km panel
restore

* (3)
reghdfe ndvi c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.seila_total temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (4)
reghdfe treecover c.(trt_overall_road trt_overall_irrigation trt_overall_else)##c.seila_total temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (5)
reghdfe ndvi trt_overall_road trt_overall_irrigation trt_overall_else c.year#c.(bombing_dummy burial_dummy memorial_dummy prison_dummy distance_to_city distance_to_road) temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (6)
reghdfe treecover trt_overall_road trt_overall_irrigation trt_overall_else c.year#c.(bombing_dummy burial_dummy memorial_dummy prison_dummy distance_to_city distance_to_road) temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (7)
reghdfe ndvi trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw = cell_count_30m], absorb(commune_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Commune, "Climate Controls", Y)

* (8)
reghdfe treecover trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw = cell_count_30m], absorb(commune_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Full, "Year FEs", Y, "Grid cell FEs", Commune, "Climate Controls", Y)

* (9)
reghdfe ndvi trt_overall_road temperature precipitation [aw = cell_count_30m] if max_trt_overall_road>0, absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Roads only, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (10)
reghdfe treecover trt_overall_road temperature precipitation [aw = cell_count_30m] if max_trt_overall_road>0, absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Roads only, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (11)
reghdfe ndvi trt_overall_irrigation temperature precipitation [aw = cell_count_30m] if max_trt_overall_irrigation>0, absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Sample", Irrigation only, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)

* (12)
reghdfe treecover trt_overall_irrigation temperature precipitation [aw = cell_count_30m] if max_trt_overall_irrigation>0, absorb(cell_id year) cluster(commune_id year)
outreg2 using table5.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Sample", Irrigation only, "Year FEs", Y, "Grid cell FEs", Y, "Climate Controls", Y)


*** Figure 6 produced in QGIS using PID dataset and commune-level administrative boundaries data


*** Figure 7

* a)
reghdfe ndvi i.cut_trt_overall_road years_since_first_road temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)

coefplot, keep(*.cut_trt_overall_road) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Num. completed roads projects") title("Treatment effects of roads on NDVI") ytitle("Effect on NDVI") ylabel(-0.02(0.02)0.08) yscale(range(-0.02(0.02)0.08)) saving(figure7_a, replace)

* b)
reghdfe ndvi i.cut_trt_overall_irrigation years_since_first_irrigation temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)

coefplot, keep(*.cut_trt_overall_irrigation) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Num. completed irrigation projects") title("Treatment effects of irrigation on NDVI") ytitle("Effect on NDVI") ylabel(-0.02(0.02)0.08) yscale(range(-0.02(0.02)0.08)) saving(figure7_b, replace)

* c)
reghdfe treecover i.cut_trt_overall_road years_since_first_road temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)

coefplot, keep(*.cut_trt_overall_road) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Num. completed roads projects") title("Treatment effects of roads on tree cover") ytitle("Effect on tree cover") ylabel(-0.1(0.1)0.3) yscale(range(-0.1(0.1)0.3)) saving(figure7_c, replace)

* d)
reghdfe treecover i.cut_trt_overall_irrigation years_since_first_irrigation temperature precipitation [aw = cell_count_30m], absorb(cell_id year) cluster(commune_id year)

coefplot, keep(*.cut_trt_overall_irrigation) vertical yline(0) graphregion(color(white)) legend(off) xtitle("Num. completed irrigation projects") title("Treatment effects of irrigation on tree cover") ytitle("Effect on tree cover") ylabel(-0.1(0.1)0.3) yscale(range(-0.1(0.1)0.3)) saving(figure7_d, replace)


********** Commune-Level Tables & Figures **********

* load the commune-level panel
use commune_panel, clear

* aggregate the distance-based road treatment measures to a single count of road projects
egen trt_overall_road = rowtotal(trt_road_?kmband)
* aggregate the distance-based irrigation treatment measures to a single count of irrigation projects
egen trt_overall_irrigation = rowtotal(trt_irrigation_?kmband)
* aggregate the distance-based treatment measures to a single count of projects, summing across project types
egen trt_overall = rowtotal(trt_?kmband)

* create an "other" project count. This is just the total project count minus the number of road and irrigation projects
gen trt_overall_else = trt_overall - (trt_overall_road + trt_overall_irrigation)

* maximum number of treatments for each cell
egen max_trt_overall = max(trt_overall), by(commune_id)
* rtopcode the max treatment measure to 20
egen cut_max_trt_overall = cut(max_trt_overall), at(0(1)20 1000) label


** Table 4

* Table 4 model (1) is run using the 30m grid cell panel, which is omitted because the dataset requires >100gb of memory to process

* (2)
reghdfe ndvi trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw=cell_count_30m], absorb(commune_id year) cluster(district_id year)
outreg2 using table4.doc, replace noni nocons drop(temperature precipitation) ctitle(NDVI) addtext("Year FEs", Y, "Spatial FEs", Commune, "Climate Controls", Y)

* (3)
reghdfe treecover trt_overall_road trt_overall_irrigation trt_overall_else temperature precipitation [aw=cell_count_30m], absorb(commune_id year) cluster(district_id year)
outreg2 using table4.doc, append noni nocons drop(temperature precipitation) ctitle(Tree cover) addtext("Year FEs", Y, "Spatial FEs", Commune, "Climate Controls", Y)


** Figure 8

* generate # of treatment labels for figure axes
forv t = 1/20 {
	loc t_lab "`t_lab' `t'.cut_max_trt_overall#c.year = `t'"
}

* a)
reghdfe ndvi ibn.cut_max_trt_overall#c.year if year<2003 & cut_max_trt_overall>0 [aw = cell_count_30m], absorb(commune_id) cluster(district_id year)
loc ref = _b[1.cut_max_trt_overall#c.year]
coefplot, keep(*.cut_max_trt_overall*) vertical yline(`ref') graphregion(color(white)) saving(figure8_a, replace) legend(off) xtitle("Total projects, 2003-2018") ytitle("Pre-2003 trend in NDVI") rename(`t_lab') xlabel(, labsize(vsmall))

* b)
reghdfe treecover ibn.cut_max_trt_overall#c.year if year<2003 & cut_max_trt_overall>0 [aw = cell_count_30m], absorb(commune_id) cluster(district_id year)
loc ref = _b[1.cut_max_trt_overall#c.year]
coefplot, keep(*.cut_max_trt_overall*) vertical yline(`ref') graphregion(color(white)) saving(figure8_b, replace) legend(off) xtitle("Total projects, 2003-2018") ytitle("Pre-2003 trend in tree cover") rename(`t_lab') xlabel(, labsize(vsmall))


********** Clean Up Directory **********

shell rm *.txt















